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Abstract 

There  has  been  growing  interest  in  the  modeling  of  proton  exchange  membrane  fuel  cells  (PEMFC)  over  the  last  two  decades.  While  a  variety  of 
steady-state  models  have  been  proposed,  literature  is  scarce  in  PEMFC  dynamic  models  and  transient  studies.  Typical  dynamic  models  for  PEM 
fuel  cell  are  empirical  current- voltage  relationships.  The  internal  transients  associated  with  reactant  and  product  species  and  other  components 
are  usually  neglected.  A  detailed  dynamic  model  for  spherical  agglomerate  in  PEM  fuel  cell  is  presented  in  this  work.  The  dynamic  model 
includes  detailed  mathematical  equations  for  conservation  of  oxygen  and  hydrogen  ions  inside  the  agglomerate.  The  agglomerate  dynamic  model 
is  simulated  for  typical  operating  conditions  inside  the  PEMFC  catalyst  layer.  Simulation  studies  show  that  the  time  scales  in  which  the  dynamics 
of  agglomerate  potential  and  concentration  of  dissolved  oxygen  respond  differ  by  several  orders  of  magnitude.  Transient  response  of  agglomerate 
current  to  step  changes  in  surface  boundary  conditions  are  also  presented.  Reasons  for  the  typical  characteristics  observed  in  the  dynamic  behavior 
of  agglomerate  current  are  also  highlighted. 

©  2005  Published  by  Elsevier  B.V. 
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1.  Introduction 

There  has  been  considerable  interest  in  the  modeling  of  pro¬ 
ton  exchange  membrane  fuel  cells  (PEMFC),  and  a  number  of 
PEMFC  models  have  been  proposed  in  the  literature  over  the 
last  two  decades.  The  seminal  works  in  PEMFC  modeling  are 
published  by  Springer  et  al.  [1]  and  Bernardi  and  Verbrugge 
[2,3].  Models  proposed  during  the  early  years  are  typically 
one-dimensional  and  accounted  for  steady-state  mass  transport 
and  electrochemical  kinetics.  Subsequently,  both  simplified  and 
complex  models  in  terms  of  dimensionality  and  physicochemi¬ 
cal  phenomena  have  been  studied.  These  models  have  been  used 
for  a  variety  of  purposes  such  as,  prediction  of  the  typical  char¬ 
acteristic  (current-potential)  curves,  parameter  and  operating 
conditions  sensitivity  analysis,  structural  and  process  optimiza¬ 
tion  studies,  and  three-dimensional  temperature,  pressure  and 
species  concentration  distributions  in  the  case  of  fuel  cell  stacks. 
A  common  feature  in  most  of  these  models  was  that  the  reac¬ 
tion  or  catalyst  layer  is  not  modeled  in  detail.  The  reaction  layer 
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is  treated  as  an  ultra-thin  layer,  thus  neglecting  the  transport  of 
reactant  gases  and  products.  Hence,  the  catalyst  layer  is  treated 
as  a  source/sink  boundary  condition  for  transport  equations  in 
the  gas  diffusion  layer.  Contrary  to  this  assumption,  even  if  gas 
phase  transport  is  neglected  on  the  consideration  of  ultra-thin 
layer,  the  presence  of  ionomer  in  the  reaction  layer  along  with 
carbon  and  platinum  makes  transport  within  the  pores  of  the 
ionomer  important.  Moreover,  catalyst  layer  is  the  region  where 
various  limiting  mechanisms  can  occur  and  thus,  can  have  a 
strong  influence  on  the  overall  performance  of  the  cell. 

It  is  widely  accepted  that  gas-diffusion  electrodes  (GDEs) 
used  in  fuel  cells  are  three-phase  electrodes  with  a  complex 
geometry  consisting  of  conducting  solid  material  for  electron 
transfer,  hydrophilic  phase  for  ion  transfer  and  open  pores  for  gas 
transport.  It  is  generally  accepted  that  the  major  difficulty  in  fuel 
cell  development  lies  with  the  kinetics  of  the  oxygen  electrode 
[4,5].  The  hydrogen  electrode  kinetics  is  relatively  much  faster 
than  that  of  the  oxygen  electrode.  The  structure  and  operation  of 
hydrophobic  electrodes  have  been  studied  for  the  first  time  under 
electron  micrographs  by  Tantram  and  Tseung  [6].  They  estab¬ 
lished  that  the  catalyst  is  present  largely  in  the  form  of  porous 
aggregates,  which  form  a  three-dimensional  network  throughout 
the  electrode.  Intermingled  with  this,  there  exists  an  interlocking 
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effective  area  of  catalyst  site  per  unit  volume  of 
agglomerate  (cm2  Pt  cm-3  agglomerate) 
activity  of  water  in  the  ionomer 
effective  area  of  Pt  per  unit  weight  of  Pt 
(cm2  Ptmg  Pt-1) 

capacitance  per  unit  volume  of  agglomerate 

(farad  cm-3  agglomerate) 

concentration  of  dissolved  oxygen  inside  the 

agglomerate  (gmol  cm-3) 

concentration  of  dissolved  oxygen  inside  the 

ionomer  film  (gmol  cm-3) 

saturation  concentration  of  oxygen  inside  the 

ionomer  pores  (gmol  cm-3) 

diffusivity  of  dissolved  oxygen  in  ionomer  pores 

inside  the  agglomerate  (cm2  s-1) 

diffusivity  of  dissolved  oxygen  in  ionomer  pores 

(cm2  s-1) 

Faraday  constant  (C  gm  equiv.-1) 

Henry’s  constant  for  air-ionomer  interface 

(atm  cm3  gmol-1) 

local  current  density  (A  cm-2  Pt) 

exchange  current  density  for  oxygen  reduction  on 

Pt  (A  cm-2  Pt) 

instantaneous  current  from  an  agglomerate  (A) 
charge  flux  inside  agglomerate  (A  cm-2) 
charge  flux  inside  ionomer  film  (A  cm-2) 
loading  of  catalyst  in  cathode  catalyst  layer 
(mg  Pt  cm-2) 

no.  of  electrons  taking  part  in  oxygen  reduction 
reaction 

partial  pressure  of  oxygen  at  the  surface  (atm) 
partial  pressure  of  water  vapor  at  the  surface  (atm) 
saturation  pressure  of  water  vapor  (atm) 

universal  gas  constant  (J  gmol  K-1) 

rate  of  oxygen  reduction  reaction  per  unit  volume 

of  agglomerate  (gmol  cm-3  s) 

radial  distance  inside  an  agglomerate  from  the 

center  (cm) 

radius  of  agglomerate  (cm) 
agglomerate  temperature  (K) 
time  (s) 

thickness  of  cathode  catalyst  layer  (cm) 
potential  of  the  solid  phase  inside  the  agglomerate 
measure  against  SHE  (V) 


Greek  symbols 

a  transfer  coefficient 

<4lm  thickness  of  ionomer  film  covering  the  agglomer¬ 
ate  (cm) 

£agg  fraction  of  volume  occupied  by  ionomer  inside 
agglomerate 

£r  fraction  of  void  volume  inside  cathode  catalyst 

layer 


Or  local  overpotential  inside  the  agglomerate  (V) 
0agg  local  potential  inside  agglomerate  (V) 

0film  local  potential  inside  the  ionomer  film  (V) 

0r  potential  at  the  ionomer  surface  (V) 

/cagg  conductivity  of  ionomer  inside  the  agglomerate 
(mho  m1 ) 

/Cmem  conductivity  of  ionomer  film  (mho  m- 1 ) 

A  water  water  content  inside  the  ionomer 

(molH20mol-1  SO^~) 


network  of  porous  hydrophobic  polytetrafluoroethylene  (PTFE). 
In  the  cell,  the  hydrophilic  catalyst  network  becomes  flooded 
with  the  electrolyte  while  the  hydrophobic  PTFE  allows  for  the 
gas  to  diffuse  through  it.  They  had  shown  that  the  performance 
of  electrodes  is  dependent  on  the  microstructure  of  porous  cat¬ 
alyst  aggregates.  However,  it  is  not  clear  whether  any  specific 
geometry  was  suggested  for  the  catalyst  aggregates. 

As  GDEs  are  difficult  to  characterize,  one  of  the  first  assump¬ 
tions  that  was  made  to  model  them  was  the  concept  of  “flooded 
agglomerates”,  introduced  by  Giner  and  Hunter  [7].  They  have 
considered  cyclindrical  geometry  for  the  agglomerates.  Results 
were  presented  for  alkaline  oxygen  electrode.  The  potential  drop 
was  assumed  to  change  only  in  the  axial  direction  and  diffusion 
of  the  dissolved  gases  was  assumed  to  be  in  the  radial  direction 
of  the  agglomerates.  The  effect  of  the  cylindrical  agglomerate 
radius  on  the  current  generated  and  its  distribution  were  studied. 
Porous  GDEs  with  the  same  assumptions  have  been  extended  to 
model  the  phosphoric  acid  fuel  cells  (PAFC)  cathode  and  anode 
in  detail  [8,5].  Various  transport  and  kinetic  processes  which 
take  place  in  the  porous  electrodes  were  taken  into  account.  The 
model  was  used  in  the  simulation  mode  for  predictive  analysis 
and  for  electrode  development  process. 

One  of  the  drawbacks  with  the  above-proposed  flooded- 
agglomerate  model  is  that  it  does  not  consider  any  tortuosity  for 
the  gas  phase  transport  as  the  agglomerates  are  assumed  to  com¬ 
pletely  extend  from  the  gas  side  to  the  electrolyte  side.  The  cylin¬ 
drical  flooded-agglomerate  model  was  modified  by  Celiker  et  al. 
[9]  considering  spherical  geometry.  They  have  investigated  their 
model  predictions  by  considering  the  cathodic  reduction  of  oxy¬ 
gen  in  alkaline  medium.  Subsequently,  many  studies  conducted 
by  various  researchers  with  the  spherical  flooded-agglomerate 
model  were  presented  for  alkaline  fuel  cells  (AFC)  [10,11]  and 
PAFC  [11-13]. 

Even  in  the  case  of  PEM  fuel  cells,  researchers  have  studied 
the  effect  of  various  phenomena  in  the  catalyst  layer  based  on 
flooded-agglomerate  model  [11,14-18].  Perry  et  al.  [11]  have 
developed  a  model  for  gas  diffusion  electrode  that  can  be  used 
as  diagnostic  tool  for  designing  of  fuel  cells.  The  model  is  one¬ 
dimensional  model  for  mass  transport  in  the  zone  where  Tafel 
kinetics  is  valid.  The  models  presented  were  generally  valid 
for  any  GDE  with  either  liquid  electrolyte  (AFC  and  PAFC) 
or  ion-exchange  membrane  (PEMFC).  The  model  was  used  to 
study  the  effects  of  mass-transport  limitations  on  the  polariza¬ 
tion  characteristics  of  oxygen  reduction  reaction  in  the  cathode. 
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Using  negligible  mass  transfer  resistance  in  the  gas  phase  the 
model  develops  a  function  for  evaluating  the  Thiele  modulus 
for  the  catalyst-binder  agglomerates.  These  relations  along  with 
ion  transport  equations  are  combined  to  develop  a  single  variable 
second-order  differential  equation.  Based  on  it,  asymptotic  solu¬ 
tions  were  developed  at  different  limiting  conditions.  The  model 
also  predicts  different  Tafel  slopes  for  distinct  regions.  The  au¬ 
thors  have  also  shown  how  the  results  may  be  used  as  a  diag¬ 
nostic  tool  for  analyzing  fuel  cell  cathode  data.  Siegel  et  al.  [14] 
have  proposed  a  steady-state  two-dimensional  PEMFC  model 
based  on  agglomerate  geometry  for  the  catalyst  layer.  The  ag¬ 
glomerates  are  characterized  by  mean  diameter  and  a  character¬ 
istic  length.  Based  on  the  model  results,  it  has  been  highlighted 
that  the  fuel  cell  performance  is  highly  dependent  on  catalyst 
structure.  Recently,  a  thin-film  agglomerate  model  with  cylin¬ 
drical  geometry  was  proposed  for  the  cathode  of  a  PEM  fuel  cell 
[15].  The  model  described  was  a  two-phase,  one-dimensional 
steady  state,  isothermal  model  for  a  PEMFC  cathode.  The  au¬ 
thors  presented  simulation  results  confirming  water- flooding  sit¬ 
uation  in  catalyst  layers.  Since  the  agglomerates  are  considered 
to  be  cylindrical  and  extend  from  gas  diffusion  layer  to  mem¬ 
brane,  the  model  assumes  that  gas  phase  transport  inside  catalyst 
layer  is  straight.  Wang  et  al.  [18]  have  investigated  transport  and 
reaction  kinetics  in  spherical  agglomerates  of  cathode  catalyst 
layer.  They  have  considered  two  types  of  spherical  agglomerates: 
the  first  one  consisting  of  a  mixture  of  carbon/catalyst  particles 
and  perfluorosulfonated  ionomer  (PFSI)  and  the  other  type  con¬ 
sisting  of  carbon/catalyst  particles  and  water- filled  pores.  The 
model  has  been  used  to  study  current  conversion,  reactant  and 
current  distribution  and  catalyst  utilization.  The  significance  of 
these  results  for  optimization  of  catalyst  layers  have  also  been 
highlighted. 

One  of  the  probable  reasons  for  neglecting  reaction  layer 
in  PEMFC  models  published  in  the  beginning  could  be  lack 
of  instrumentation  to  characterize  its  morphology  accurately. 
With  the  availability  of  advanced  microscopy  instruments  like 
scanning  electron  microscopy  (SEM)  and  transmission  electron 
microscopy  (TEM),  researchers  have  been  able  to  study  the 
morphology  of  complex  nanostructures  such  as,  PEM  fuel 
cell  electrodes.  Middleman  [19]  has  studied  the  structure  of 
membrane-electrode-assembly  (MEA)  using  high-resolution 


scanning  electron  microscopy  (HR-SEM).  In  his  investigation 
he  has  shown  that  the  catalyst  layer  consists  of  a  random  distri¬ 
bution  of  pores  and  particles,  see  Fig.  1 .  It  was  also  shown,  using 
higher  magnification,  that  the  particles  are  agglomerates  of  much 
smaller  particles  coated  with  a  film  of  Nafion.  From  the  images  in 
the  figure  it  can  be  clearly  seen  that  the  agglomerates  are  spheri¬ 
cal  in  shape.  Lee  et  al.  [20]  and  Liu  et  al.  [21]  have  also  published 
their  investigations  of  PEM  fuel  cell  electrodes  using  SEM/TEM 
that  corroborate  Middleman’s  work.  Therefore,  spherical  ag¬ 
glomerates  are  not  an  abstract  conceptualization,  but  a  realistic 
representation  of  cathode  catalyst  layers  in  PEM  fuel  cells  for 
modeling  purposes. 

In  this  work,  the  transient  characteristics  of  various  trans¬ 
port  and  electrochemical  phenomena  are  studied  in  spheri¬ 
cal  agglomerates  using  its  dynamic  model.  For  PEM  fuel 
cells,  flooded-agglomerate  is  a  uniform  mixture  of  Pt  sup¬ 
ported  carbon  particles  and  ionomer.  In  addition,  to  the  study 
the  effect  of  ionomer  layer,  we  also  consider  a  thin  film  of 
ionomer  covering  the  agglomerate.  We  propose  a  detailed  dy¬ 
namic  study  of  the  spherical  agglomerate  for  the  following 
reasons: 

•  This  is  the  first  time  in  PEM  fuel  cell  literature  that  an  agglom¬ 
erate  dynamic  model  is  being  studied.  The  dynamic  model  has 
been  used  to  study  the  dynamics  of  both  dissolved  concentra¬ 
tion  of  oxygen  (co2,agg  and  co2,film)  and  hydrogen  ions  (0agg 
and  0mm). 

•  An  agglomerate  is  a  microcosm  of  the  reaction  layer.  There 
are  species  transport  and  electrochemical  phenomena  taking 
place  inside  the  agglomerate.  Transport  is  in  the  form  of  dif¬ 
fusion  of  dissolved  oxygen  and  hydrogen  ions  through  the 
water-filled  ionomer  pores,  and  oxygen  reduction  reaction 
at  the  active  catalyst  sites  producing  water.  Hence,  it  is  ex¬ 
pected  that  the  steady-state  and  dynamic  characteristics  of  an 
agglomerate  are  similar  to  that  of  the  fuel  cell  characteristics. 

•  The  dynamics  of  a  complete  PEM  fuel  are  not  very  well  un¬ 
derstood.  Moreover,  the  final  dynamic  model  for  the  electrode 
is  likely  to  be  complicated.  In  view  of  this,  we  intend  to  study 
the  simplifications  that  are  possible  for  the  complete  reac¬ 
tion  layer  level  model  using  the  proposed  dynamic  spherical 
agglomerate  model. 


Higher-magnification  SEM  image,  showing 
that  the  particles  are  in  fact  agglomerates  of 
much  smaller  particles 


Still  higher-magnification  SEM  image, 
showing  that  the  particles  are  in  agglomerates 
of  much  smaller  particles 


Higher  electron  energy  in  the 
HR-SEM  shows  that  the  agglomerates 
are  built  up  from  30  nm  carbon 
particles,  and  that  these  particles  are 
coated  with  finely  dispersed  platinum 


Fig.  1.  Scanning  electron  microscope  (SEM)  images  of  PEM  fuel  cell  electrode  [19]. 
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•  We  will  demonstrate  very  interesting  dynamic  behavior  such 
as  phenomenon  occurring  at  different  time  scales  and  delays 
at  the  agglomerate  level.  The  insights  developed  through  this 
model  will  significantly  enhance  the  analysis  of  the  full  dy¬ 
namic  electrode  model. 

•  The  proposed  dynamic  spherical  agglomerate  model  is  the 
back-bone  for  an  integrated  reaction  layer  model  for  the  cath¬ 
ode  and  anode. 

The  dynamic  model  for  the  agglomerate  describes  the  conserva¬ 
tion  of  dissolved  oxygen  and  hydrogen  ions,  which  is  described 
in  detail  in  the  subsequent  section.  In  the  subsequent  sections 
of  the  paper,  hydrogen  ions  and  charge  are  used  synonymously. 
Typical  steady- state  agglomerate  i-v  characteristic  curves  are 
studied.  Simulation  results  to  show  the  effect  of  boundary  con¬ 
ditions  and  model  parameters  on  the  i-v  curves  are  also  pre¬ 
sented.  The  dynamic  model  is  simulated  to  show  the  effect  of 
step  changes  in  the  surface  boundary  conditions  on  the  agglom¬ 
erate  current.  Characteristic  behavior  observed  through  simula¬ 
tions  in  the  agglomerate  current  are  also  explained. 


Carbon 


Thin  film  of 
ionomer  (5) 


Platinum 
Nanoparticles 


Hydrated 

Ionomer 


Fig.  2.  Spherical  flooded-agglomerate  with  ionomer  film. 


•  potential  of  the  solid  phase  (C  +  Pt)  is  assumed  to  be  uniform 
throughout  the  agglomerate. 


2.  Spherical  agglomerate  dynamic  model 

A  schematic  of  the  spherical  agglomerate  considered  in  this 
work  is  illustrated  in  Fig.  2.  The  agglomerate  consists  of  a  thin 
film  of  ionomer  covering  a  uniform  mixture  of  carbon  supported 
Pt  nanoparticles  and  ionomer  in  between  the  carbon  agglomer¬ 
ates.  At  the  surface,  oxygen  present  in  the  gas  dissolves  into 
the  water  present  inside  the  pores  of  the  ionomer.  At  the  sur¬ 
face  it  is  assumed  that  there  exists  an  equilibrium  between  the 
partial  pressure  of  oxygen  in  the  gas  phase  and  the  dissolved  con¬ 
centration  in  the  ionomer  phase.  The  dissolved  oxygen  diffuses 
through  the  ionomer  pores  towards  the  center  of  the  agglom¬ 
erate.  No  reactions  are  considered  in  the  thin  film  of  ionomer. 
Within  the  uniform  mixture  region,  dissolved  oxygen  is  assumed 
to  diffuse  towards  the  center  simultaneously  undergoing  the  fol¬ 
lowing  chemical  reaction  at  the  active  catalyst  sites: 

02  +  4H+  +  4e~  -*  2H20  (1) 


The  conservation  equations  for  oxygen  and  hydrogen  ions  are 
written  separately  for  the  ionomer  film  and  the  agglomerate. 


2.7.  Conservation  equations  inside  ionomer  film 


Since  no  reactions  are  considered  inside  the  ionomer  film, 
the  dynamic  conservation  equation  for  oxygen  inside  the  film 
can  be  written  as 


d 

dt 


(c02,film) 


mem 


d 


/  2  dCQ2,film\ 
V  dr  ) 


where  r  is  the  radial  distance  from  the  center  of  the  agglomer¬ 
ate  and  7)o2,mem  is  the  diffusivity  of  dissolved  oxygen  inside 
the  ionomer  pores.  With  the  above  assumptions,  the  dynamic 
conservation  equation  for  charge  flux  would  reduce  to 


-V  •  jfilm  =  0 


The  hydrogen  ions  produced  in  the  anode  catalyst  layer  travel 
across  the  membrane  layer  and  reach  active  catalyst  sites  inside 
the  agglomerate  through  a  network  of  micropores  in  the  ionomer. 
The  following  assumptions  are  considered  for  setting  up  the 
dynamic  model  equations: 


where  the  charge  flux,  jmm,  is  related  by  ohm’s  law  to  the  po¬ 
tential  inside  film  (0mm)  as 

jfilm  —  ^inciri  ^  (0film)  (4) 

Therefore,  Eq.  (3)  can  be  written  as 


•  axisymmetric  conditions  are  assumed  for  concentration  and 
charge  distribution  inside  the  agglomerate, 

•  pressure  and  temperature  conditions  are  uniform  throughout 
the  agglomerate, 

•  the  contribution  due  to  convection  for  species  flux  is  negligi¬ 
ble, 

•  charge  accumulation  takes  place  only  near  the  surface  of  Pt, 

•  Butler- Volmer  kinetics  are  considered  for  the  oxygen  reduc¬ 
tion  reaction, 

•  physical  properties  of  the  ionomer  are  considered  same  as 
that  of  the  membrane, 


^incm  9 
r2 


The  membrane  conductivity,  /cmem,  is  a  function  of  water  content 
inside  the  ionomer  pores,  Awater,  which  in  turn  is  a  function  of 
the  activity  of  water.  The  activity  is  water  is  given  by 


<2h2o  = 


PH2Q,r 

Ph2o 


where  the  saturation  pressure  of  water  vapor,  ph2Q  *s  depen¬ 
dent  on  the  local  temperature  conditions  and  is  calculated  by  the 
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following  equation: 

pfc,0  =  ^(7'.ggr4’28310l23MI,“c’37-8/r*«))  (7) 

The  ionomer  water  content  is  calculated  from  [1]: 

7-water  —  0.043  +  17.81<2h2o  —  39.85 ciii2Q  +  36.0<2^j2q, 
forO  <  aii2o  <  1 

=  14.0  +  1.4(<2h2o  —  1), 

for  1  <  au2o  <  3  (8) 


The  membrane  conductivity  is  calculated  by  the  following  equa¬ 
tion  [1]: 


^mem  —  (0.005 139A water  -  0.00326) 

1  1 


x  exp 


1268.0 


303.0  Z 


agg 


(9) 


2.2.  Conservation  equations  inside  the  agglomerate  region 


The  dynamic  conservation  equation  for  oxygen  inside  the 
agglomerate  can  be  written  as 


^  \  _  ^°2, agg 

02,aggJ-  r2 


3co2, 


agg 


dr 


+  R 


orr 


(10) 


where  Do2,agg  is  the  diffusion  coefficient  of  dissolved  oxygen 
in  the  ionomer  pores  inside  the  agglomerate.  It  is  related  to 
Ao2,mem  and  fraction  of  ionomer  volume  inside  the  agglom¬ 
erate  (£agg)  as 

^02,agg  —  £agg^02,mem  (H) 

The  oxygen  reduction  reaction  ( Ron )  is  given  by  Butler- Volmer 
kinetics  as 


a  friz 
nF 


(12a) 


aa 


lo 


c02,agg 


C 


exp 


o2 


(  anQqvF 
\  R1 agg 


—  exp 


f  (1  ~  a)nQqYF\  1 

V  RT&gg  J  J 


(12b) 


The  local  overpotential  (qr)  appearing  in  the  Butler- Volmer  ki¬ 
netics  is  the  defined  by  the  following  equation 


Or  =  A0(s,  agg)  -  A0e(s,  agg)  (13a) 

—  {0s  —  0agg}  —  {0e,s  —  0e,agg}  (13b) 

where  0S  is  the  potential  of  the  solid  phase  and  0agg  is  ionomer 
potential  adjacent  to  carbon  and  platinum.  The  subscript  e  de¬ 
notes  equilibrium  conditions.  If  oxygen  electrode  is  treated  as 
the  reference  electrode,  then  solid  phase  potential  0S  =  0  and 
A0e(s,  agg)  =  0.  In  such  a  case,  the  local  overpotential  appear¬ 
ing  in  the  Butler- Volmer  equation  (qY)  is  numerically  equal  to 
the  local  ionomer  potential  (0agg).  The  active  surface  area  of 
platinum  per  unit  volume  of  agglomerate  (aa)  is  related  to  plat¬ 
inum  loading  (rapt),  effective  area  of  platinum  (apt),  catalyst 


layer  thickness  (£rl),  and  void  fraction  inside  the  catalyst  layer 
(£r)  by  the  following  equation 


^pt^pt 
(1  —  £r)tRL 


(14) 


The  saturation  concentration  of  dissolved  oxygen,  Cq2  ,  is  given 
by  Henry’s  law  applied  at  the  surface  of  the  ionomer  film: 


77 -  P02,r  (15) 

^02,mem  / 

Writing  the  dynamic  conservation  equation  for  charge  gives 
d 

—  (CaggOr)  —  —  V  ‘jagg  T  HqFR0YY  (16) 

ot 

where  qY  is  the  local  overpotential  inside  the  agglomerate  at  a 
radial  distance  r  from  the  center  and  Cagg  is  capacitance  per 
unit  volume  of  agglomerate.  For  writing  the  above  conserva¬ 
tion  equation,  both  capacitative  and  faradaic  currents  are  taken 
into  consideration.  The  left-hand  side  of  the  equation  represents 
accumulation  of  hydrogen  ions  in  the  double  layer  near  the  plat¬ 
inum  surface.  The  faradaic  current  is  represented  by  the  reaction 
term  (R0rr)-  The  charge  flux  inside  the  agglomerate,  jagg,  is  given 
by 


jagg  —  ^aggV(0agg) 


(17) 


The  conductivity  of  hydrated  ionomer  inside  the  agglomerate, 
/cagg,  is  related  to  membrane  conductivity  as 


k  -  p2/2k 
^agg  —  fcaggKmem 

Eq.  (16)  can  be  written  as 


(18) 


d 

Jt 


((-agg  Or) 


/Cagg  ^  (2  90agg\ 

r2  dr  \  dr  J 


T  RqFRqt^ 


(19) 


Eqs.  (2),  (5),  (10)  and  (19)  are  solved  for  the  concentration  of 
dissolved  O2  and  potential  inside  the  agglomerate  and  ionomer 
film.  The  initial  and  boundary  conditions  are  given  in  Table  1 . 
Henceforth,  the  potential  at  the  surface  of  ionomer  film  (0r)  will 
be  called  as  surface  potential  in  the  paper. 

The  agglomerate  current  can  be  calculated  by  several  meth¬ 
ods.  For  the  steady-state  model,  the  agglomerate  current  was 
calculated  by  the  following  equations,  which  are  all  equivalent 


j  j-t  a  2  r\  9co2,film 

7agg  —  ^eE47rragg/9o2,film  lf=ragg 

(20a) 

—  n  F4irr2  Dn  ^c°2,agg , 

“  nQr^7zr^ggu  02,agg  ^  lr=ragg 

(20b) 

A  2  ^0film  | 

-  47Tragg/':mem  ^  1  r=ragg 

(20c) 

-  4nr2  k  9<^aggl 

-  ^^agg^agg  ^  lr=ragg 

(20d) 

fraSB  0 

=  /  47tr2(nQFROYY)dr 

Jo 

(20e) 

In  the  case  of  dynamic  model,  for  determining  the  instanta¬ 
neous  current  generated  inside  an  agglomerate,  local  reaction 
rate  along  with  charge  accumulation  term  is  integrated  over  the 

R.  Madhusudana  Rao,  R.  Rengaswamy  /  Journal  of  Power  Sources  158  (2006)  110-123 


Table  1 

Initial  and  boundary  conditions  for  spherical  agglomerate  in  Fig.  2 
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Number 


Initial  and  boundary  conditions 


t  =  0 
IC1 
IC2 
IC3 
IC4 

r  =  ragg  T  ^film 
BC1 


C02,filmV/' 

C02,aggVr 

0filmVr 

0aggV/* 


C02,film 


1 


H0 


2,  mem 


P02,rVt 


BC2 

0film  —  0r  Vf 

=  ragg 

BC3 

c02,agg  =  £  agg  C02,  film  W 

BC4 

n  9co2,film 

*A>2 .  mem  ^  —  ^02,agg 

BC5 

0agg  —  0filir A/ 

BC6 

90film  90agg 

/Cmem  9r  “  ^  9r  W 

=  0 

9c°2.agg  =()  yt 
dr 

BC7 

BC8 

9^gg  =  0  Vf 

dr 

dc02, 


agg 


dr 


Vt 


Comments 


Obtained  from  steady-state  solution 
Obtained  from  steady-state  solution 
Obtained  from  steady-state  solution 
Obtained  from  steady-state  solution 


Equilibrium  at  interface 
Surface  boundary  condition 

Continuity 

Species  flux  condition 
Continuity 

Charge  flux  condition 


entire  agglomerate  volume  and  expressed  in  current  equivalent 
and  is  given  by  the  following  equation: 

Tagg  ?  f  9  I 

^agg  —  J  4-Jir  ^ —  HqFRq^  +  ~(^r^agg)  j  dr  (21) 

3.  Results  and  discussion 

The  model  equations  were  set  up  Maple  (version  9)  and  were 
solved  using  fsolve  (for  steady-state  case)  and  ode  15s  (for  dy¬ 
namic  case),  which  are  available  in  MATLAB.  For  generating 
initial  guess  for  the  dynamic  model,  steady-state  model  corre¬ 
sponding  to  the  dynamic  model  equations  enumerated  above  are 
solved  first. 

3.1.  Steady-state  model  simulation 

The  steady-state  model  is  set  up  such  that,  surface  potential 
(0r),  partial  pressures  of  oxygen  (/?o2,r)  and  water  vapor  (pn2o,r) 
at  the  surface  of  ionomer  film  and  temperature  of  the  agglom¬ 
erate  (Tagg)  are  taken  as  inputs.  The  dissolved  concentration  of 
oxygen  (co2,agg  and  co2,film),  potential  variation  inside  ionomer 
film  and  agglomerate  (0a gg  and  0fnm),  and  the  current  generated 
from  the  agglomerate  (/agg)  are  the  output  variables.  Simulation 
studies  are  carried  out  using  the  steady- state  model  to  study  the 
i-v  characteristic  curve  and  the  influence  of  boundary  condi¬ 
tions  and  other  parameters  on  the  curve.  The  base  case  model 
parameters  and  their  range  for  parameter  sensitivity  studies  are 
listed  in  Table  2. 

For  discretizing  the  PDEs,  a  certain  number  of  points  are 
chosen  along  the  radius  of  the  agglomerate  and  finite  differ¬ 
ence  techniques  are  applied.  Before  describing  the  model  results, 
the  sensitivity  of  the  results  to  the  number  of  points  chosen  is 
briefly  discussed  here.  Preliminary  steady-state  simulation  re¬ 


sults  showed  that  oxygen  concentration  profiles  inside  ionomer 
film  and  agglomerate  were  sensitive  to  the  number  of  points 
chosen  along  the  radius.  As  a  result,  the  i-v  curves  also  showed 
variation.  Hence,  simulations  were  performed  by  increasing  the 
number  of  points  till  the  results  showed  no  variation.  In  addition, 
it  was  also  observed  that  the  number  of  points  inside  the  agglom¬ 
erate  (Nagg)  and  ionomer  film  (A^fiim)  have  to  be  same.  In  all  the 
results  presented  below,  Aagg  =  Nfi\m  =  500.  For  the  sake  of 
clarity,  even  though  500  points  are  chosen  along  the  radius  in 
the  agglomerate  and  ionomer  film,  only  fewer  number  of  points 
are  shown  in  the  plots.  Figs.  3  and  4  illustrate  the  concentration 
profiles  of  dissolved  oxygen  inside  the  agglomerate  and  ionomer 
film,  respectively,  at  different  surface  potentials.  The  plots  show 
that  as  surface  potential  is  decreased,  gradients  of  the  concen¬ 
tration  profiles  become  steeper.  This  is  due  to  an  increase  in 
reaction  rate  as  the  surface  potential  is  decreased.  When  the  sur¬ 
face  potential  is  decreased  beyond  —0.20  V,  it  was  observed  that 
the  concentration  profiles  have  more  steeper  gradients.  These 
results  show  that  even  though  the  thickness  of  ionomer  film  is 
very  small  (in  the  range  of  0.01-0.10  pm)  considerable  con¬ 
centration  gradients  exist,  contrary  to  the  assumptions  made  by 
other  researchers  [15,18]. 

Since  an  agglomerate  is  a  microcosm  of  the  reaction  layer, 
it  will  be  expected  that  the  characteristic  curve  of  an  agglomer¬ 
ate  is  similar  to  that  of  the  fuel  cell  characteristic  curve.  Fig.  5 
illustrates  the  i-v  characteristic  curve  for  the  spherical  agglom¬ 
erate,  which  conforms  to  the  typical  fuel  cell  characteristic  plot. 
The  curve  is  marked  by  the  legend  “With  variation  in  potential” 
in  Fig.  5.  The  other  i-v  curve  in  the  figure  corresponds  to  the 
case  when  potential  is  considered  constant  inside  the  ionomer 
film  and  agglomerate,  an  assumption  considered  by  Wang  et  al. 
[  1 8]  for  PFSI  filled  spherical  agglomerate.  With  this  assumption, 
only  the  steady- state  oxygen  conservation  equation  needs  to  be 
solved.  In  fact,  an  analytical  solution  can  be  obtained  for  the 
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Table  2 

Model  inputs  and  parameters 


Parameter 


Inputs 


Base  case  value 


P02J 

0.50 

PH20,r 

0.24 

0r 

0.0  to  —0.40 

Agg 

353.15 

odel  parameters 

Apt 

1000 

Urxn 

D  02,mem 

20.0  X  10“6 

3.1  x  10-3  exp 

f  2768 

\  Agg 

F 

,  mem 

96485.3 

1.33  x  106  exp  | 

(  666.0 

\  Agg 

io 

6.0  X  10“8 

rapt 

0.40 

nc 

4 

R 

8.314 

ragg 

1.0  X  1CT4 

Agg 

353.15 

tRL 

30 

a 

0.50 

^film 

1.0  x  10“5 

£agg 

0.40 

0.40 

Range  for  parameter  studies 


Units 


Reference 


0.05-0.50 

atm 

0.10-0.60 

atm 

V 

K 

_ 

cm2  Pt  mg  - 1  Pt 

[15] 

Farad  cm-2  Pt 

[22,23] 

2  -l 
cm  s 

[24,2] 

Cgequiv.-1 

atm  cm3  gmol- 1 

[3] 

A  cm-2 

[14] 

0.10-1.0 

mg  Pt  cm-2 

[15] 

J  gmol  K- 1 

1.0  X  10"5  to  1.0  X  10“4 

cm 

K 

[Jim 

5.0  x  10-6  to  5.0  x  10"5 

cm 

0.10-0.50 

[3] 

[3] 


agglomerate  current  (/agg),  which  was  obtained  using  Maple.  A 
comparison  of  the  two  i-v  curves  clearly  shows  the  difference 
and  the  significance  of  considering  variation  in  potential  inside 
the  ionomer  film  and  agglomerate  instead  of  assuming  it  to  be 
constant. 

Another  important  feature  to  be  noted  from  the  plot  is  that 
the  current  from  the  agglomerate  stagnates  beyond  a  certain 
surface  potential  (—0.35  V).  This  is  due  to  the  fact  that  at  these 
potentials,  the  dissolved  oxygen  is  completely  consumed  within 
a  short  radial  distance  from  the  surface  inside  the  agglomerate 
region.  In  fact,  it  can  be  observed  from  Fig.  3  that  beyond  a 


Fig.  3.  Normalised  concentration  of  O2  as  a  function  of  radial  distance  inside 
agglomerate  for  different  surface  potentials  (0r). 


certain  potential,  oxygen  is  completely  utilised  at  the  surface  of 
the  agglomerate  itself  and  there  is  not  enough  supply  of  oxygen 
to  reach  the  center  of  the  agglomerate.  Hence,  a  large  part  of  the 
agglomerate  does  not  produce  any  current  due  to  lack  of  oxygen. 

Parametric  studies  were  carried  out  on  the  steady- state  model 
of  the  agglomerate  to  study  the  effect  of  boundary  conditions 
and  model  parameters  on  the  characteristic  curve.  For  paramet¬ 
ric  studies  using  the  model,  a  range  of  values  are  considered 
for  a  particular  parameter  and  base  case  values  are  taken  for  the 
other  parameters.  The  steady  state  model  was  simulated  for  dif¬ 
ferent  values  of  surface  boundary  conditions  (po2j,  anc^  PH2o,r) 


Fig.  4.  Normalised  concentration  of  02  as  a  function  of  radial  distance  inside 
ionomer  film  for  different  surface  potentials  (0r). 
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Fig.  5.  I-V  characteristic  curve  for  agglomerate  shown  in  Fig.  2. 


and  model  parameters  (mpt,  £agg,  ragg  and  Fig.  6(a)  shows 
the  characteristic  curves  for  different  values  of  partial  pressure 
of  oxygen  at  the  ionomer  surface.  With  an  increase  in  partial 
pressure  of  oxygen  in  the  gas  phase,  the  equilibrium  concen¬ 
tration  of  dissolved  oxygen  just  inside  the  ionomer  increases. 
Hence,  as  it  would  be  expected,  larger  amounts  of  current  are 
generated  from  the  agglomerate  and  this  trend  can  be  clearly 
seen  in  the  figure.  Contrary  to  the  above  result  and  expectations, 
increasing  the  partial  pressure  of  water  vapor  at  the  surface  does 
not  have  any  significant  effect  on  the  i-v  characteristics  of  the 
agglomerate,  as  shown  in  Fig.  6(b).  This  can  be  explained  as 
follows:  with  an  increase  in  partial  pressure  of  water  vapor,  the 
amount  of  water  content  in  the  ionomer  increases  and  thus,  the 
ionomer  conductivity  also  increases.  This  is  expected  to  result 
in  an  increase  in  the  current,  since  current  is  given  by 

,  —  A-rrr 2  ^agg  I  _  .  2  #  mem  . 

fagg  =  -K  agg^agg  — —  lr=ragg  =  -Kmem^agg  —  I  r=ragg 

(22) 


On  the  other  hand,  when  the  conductivity  increases,  the  poten¬ 
tial  drop  from  the  surface  to  the  agglomerate-ionomer  interface 
is  less  and  thus,  the  potential  gradient  is  decreased.  These  two 
opposing  effects  nullify  each  of  their  contributions  and  this  re¬ 
sults  in  producing  no  change  in  the  agglomerate  current.  Use  of 
the  above  equation  for  calculating  agglomerate  current  does  not 
mean  that  conductivity  of  membrane  is  the  limiting  factor.  The 
current  also  depends  on  the  potential  gradient  appearing  in  the 
equation.  Since  it  is  easy  to  relate  the  effect  of  partial  pressure  of 
water  vapor  on  conductivity  of  membrane,  the  above  equation 
was  chosen  to  describe  the  effect  of  partial  pressure  of  water 
vapor  on  agglomerate  current.  Moreover,  if  we  use  the  equation 
based  on  oxygen  reaction  rate  for  current  calculation: 


4gg  — 


agg 


47Tr2(nQFRorr)  dr 


ro 


r 4mi  ^ 


exp 


an&rjvF 


RT 


agg 


exp 


(1  -  a)nGrjTF 


RT 


dr 


(23) 


agg 


the  effect  of  partial  pressure  of  water  vapor  cannot  be  easily 
explained  as  it  can  be  done  using  the  earlier  equation. 

One  of  the  limitations  of  this  result  is  that,  when  carrying  out 
the  parametric  studies,  partial  pressures  of  water  vapor  and  oxy¬ 
gen  are  treated  independently.  But  in  actual  fuel  cell  operation, 
this  will  not  the  case.  This  result  will  be  further  investigated  in 
our  research  work  on  the  dynamic  model  for  PEM  fuel  cell  cath¬ 
ode  that  is  currently  under  progress.  It  is  necessary  to  highlight 
at  this  juncture,  that  though  all  the  four  plots  in  Fig.  6(b)  seem  to 
merge,  upon  magnification  it  was  observed  that  this  was  not  the 
case.  At  a  constant  potential,  the  agglomerate  currents  differed 
from  each  other  in  third  and  fourth  decimal  points.  This  is  also 
corroborated  by  the  corresponding  dynamic  plots  (see  Fig.  16). 

The  model  was  also  simulated  for  different  agglomerate 
radii  (ragg).  The  i-v  characteristic  plots  are  shown  in  Fig.  7(a). 
At  a  constant  surface  potential,  with  increasing  agglomerate 
radii,  larger  surface  area  is  available  for  reaction.  Hence,  the 


Fig.  6.  Parametric  studies  of  agglomerate  characteristic  curve:  (a)  variation  in  po2tT  and  (b)  variation  in  pn2o,r- 
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Fig.  7.  Parametric  studies  of  agglomerate  characteristic  curve:  (a)  variation  in  ragg  and  (b)  variation  in  mpt. 


agglomerate  current  increases.  The  effect  of  change  in  platinum 
loading  (rapt)  on  i-v  characteristic  plot  is  shown  in  Fig.  7(b). 
This  figure  shows  some  unique  characteristics,  different  from 
Figs.  6(a)  and  7(a).  It  can  be  observed  that  at  a  given  potential 
(0r),  agglomerate  current  increases  with  increasing  platinum 
loading.  This  effect  in  the  agglomerate  is  clearly  seen  for 
surface  potentials  in  the  range  —0.15  to  —0.35  V.  At  surface 
potentials  between  0.0  and  —0.15  V,  the  effect  is  not  clearly 
observed  from  the  plots  due  to  the  scale  on  the  v-axis.  But  upon 
zooming  on  the  plots  between  0.0  and  —0.15  V,  the  effect  could 
be  clearly  seen.  At  surface  potentials  <  —0.35  V,  all  the  plots 
converge  to  the  same  limiting  current  value.  This  is  due  to  the 
fact  that  the  agglomerate  current  is  limited  by  the  total  amount 
of  dissolved  oxygen  and  thus,  platinum  loading  does  not  have 
any  effect.  This  result  suggests  that  there  is  an  optimum  loading 
for  the  platinum  catalyst  inside  the  PEMFC  reaction  layer. 

The  effect  of  different  ionomer  film  thickness  on  the  i-v 
curves  is  shown  in  Fig.  8(a).  At  a  particular  surface  poten¬ 
tial,  with  an  increase  in  film  thickness,  the  agglomerate  cur¬ 
rent  decreases.  This  is  due  to  the  fact  that  an  ionomer  film  with 
higher  thickness  provides  a  higher  resistance  for  diffusion  of 
dissolved  oxygen.  This  results  in  low  concentrations  of  oxygen 


at  the  agglomerate-ionomer  interface  and  hence,  less  amount 
of  reaction  leading  to  a  decrease  in  current  generated.  On  the 
other  hand,  at  a  constant  surface  potential,  when  the  fraction  of 
volume  occupied  by  ionomer  inside  active  agglomerate  region 
(€agg)  is  increased,  increase  in  agglomerate  current  is  observed, 
Fig.  8(b).  This  is  a  result  of  the  boundary  condition,  BC3,  in 
Table  1 .  Presence  of  more  ionomer  in  between  the  carbon  ag¬ 
glomerates  results  in  a  higher  supply  of  dissolved  oxygen  and 
thus,  more  reaction  and  more  current.  It  should  be  noted  that  for 
this  parameter  study,  catalyst  loading  (rapt)  and  active  catalyst 
area  (apt)  are  kept  constant  at  their  base  case  values. 

3.2.  Dynamic  model  simulation 

The  agglomerate  dynamic  model  was  simulated  for  various 
conditions  to  gain  insight  into  the  dynamic  characteristics  of  an 
agglomerate.  The  dynamic  model  is  solved  using  ode  15s,  an  ode 
solver  in  MATLAB.  MATLAB’s  Simulink  environment  is  used 
to  simulate  the  dynamic  model.  Inputs  to  the  model  and  model 
parameters  are  passed  from  the  Simulink  environment.  Fig.  9 
illustrates  the  set  up  of  the  dynamic  model  for  spherical  ag¬ 
glomerate  in  MATLAB’s  Simulink  environment.  The  model  is 


Fig.  8.  Parametric  studies  of  agglomerate  characteristic  curve:  (a)  variation  in  <5nim  and  (b)  variation  in  £a gg- 
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Fig.  9.  Agglomerate  dynamic  model  in  Matlab’s  Simulink  environment. 


set  up  such  that  ionomer  film  surface  potential  (0r),  partial  pres¬ 
sures  of  oxygen  (po2j)  and  water  vapor  (/?H?o,r)  at  the  surface 
of  ionomer  film  and  temperature  of  the  agglomerate  (Tagg)  are 
taken  as  inputs.  The  dissolved  concentration  of  oxygen  (co2,agg 
and  co2,fiim)>  potential  (</>agg  and  0mm),  and  the  current  gener¬ 
ated  from  the  agglomerate  (7agg)  are  the  output  variables.  Ini¬ 
tial  conditions  for  simulation  are  provided  by  the  corresponding 
steady- state  model. 

To  determine  the  minimum  time  for  which  the  dynamic  sim¬ 
ulation  has  to  be  run,  the  dynamic  model  equations  are  non- 
dimensionalised.  Various  time  constants  are  identified  and  it  was 
found  that  diffusive  term  of  the  potential  equation  (Eq.  (19))  has 
the  least  time  constant  for  the  given  set  of  parameters.  The  dy¬ 
namic  characteristics  of  an  agglomerate  are  studied  by  introduc¬ 
ing  a  step  in  surface  potential  and  monitoring  the  transients  in 


Fig.  10.  Dynamic  simulation  for  agglomerate  potential,  fagg'  (a)  simulation 


the  output  variables  (co2,agg,  co2,mm,  </>agg,  0fiim,  and  7agg).  Fig. 
10  shows  dynamic  simulation  results  for  variation  in  agglomer¬ 
ate  potential  (0agg)  for  a  step  in  surface  potential  from  —0.05  to 
—0.10  V.  Different  lines  in  the  figure  correspond  to  the  different 
radial  locations  inside  the  spherical  agglomerate.  Even  though 
the  agglomerate  is  divided  into  500  concentric  spherical  shells, 
profiles  for  only  1 1  radial  locations,  which  includes  r  =  0  and 
r  =  ragg,  are  illustrated  in  the  figures  for  clarity.  Fig.  10(a)  shows 
simulation  results  for  time  between  t  =  0  and  1 .4  x  10“6  s.  It 
can  be  clearly  seen  that  agglomerate  potential  reaches  steady- 
state  well  within  the  simulation  time.  On  the  other  hand,  when 
the  model  is  simulated  from  t  =  0  to  0.035  s  (Fig.  10(b)),  no 
transients  are  observed  for  all  the  radial  locations.  All  the  1 1 
lines  plotted  in  the  figure  coincide  and  an  instantaneous  step 
response  is  observed. 
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All  the  eleven  lines  are  coinciding  in  this  plot 
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Fig.  11.  Dynamic  simulation  for  potential  inside  ionomer  film,  (a)  simulation  from  time  t  =  0  to  1.4  x  10  6  s  and  (b)  simulation  from  t  =  0  to  0.035  s. 


The  corresponding  dynamic  simulation  results  for  poten¬ 
tial  inside  the  ionomer  film  (0mm)  are  shown  in  Fig.  11.  Even 
though  while  writing  the  dynamic  model,  no  charge  accumula¬ 
tion  was  considered  inside  the  ionomer  film,  it  still  exhibits  a 
small  amount  of  transients  when  the  model  is  simulated  for  the 
time  between  t  =  0  and  1.4  x  10~6  s  (Fig.  11(a)).  These  tran¬ 
sients  arise  due  to  the  boundary  conditions  at  the  agglomerate- 
ionomer  interface.  Similar  to  the  dynamic  behavior  of  0agg,  when 
the  model  is  simulated  from  t  =  0  to  0.035  s  (Fig.  11(b)),  no 
transients  are  observed. 

The  corresponding  dynamic  simulation  results  for  dissolved 
concentration  of  oxygen  inside  agglomerate  and  ionomer  film 
are  shown  in  Figs.  12  and  13,  respectively.  The  concentrations 
are  plotted  in  normalized  scale.  The  model  is  run  under  steady 
conditions  from  t  =  0  to  0.01  s  and  then  a  step  in  surface  po¬ 
tential  from  —0.05  to  — 0.10  V  is  introduced.  Fig.  12(a)  shows 
the  simulation  results  for  time  between  t  =  0  and  1.4  x  10-6  s. 
A  flat  profile  is  observed  at  all  radial  locations  for  the  complete 
time  interval.  This  can  be  explained  by  the  fact  that  the  time  in¬ 
terval  being  very  small,  the  concentrations  at  all  radial  locations 
do  not  show  any  change  from  their  steady-state  profile  corre¬ 
sponding  to  the  surface  potential  of  —0.05  V.  The  transients  in 
oxygen  concentration  are  observed  only  when  the  model  is  sim¬ 


ulated  from  t  =  0  to  0.035  s,  as  illustrated  in  Fig.  12(b).  The 
transients  exhibited  by  the  dissolved  concentration  of  oxygen 
inside  the  ionomer  film  (Fig.  13)  are  similar  to  that  of  the  con¬ 
centration  inside  the  agglomerate.  It  should  be  noted  that  the 
range  of  y-axis  in  Figs.  12  and  13  is  very  small.  This  is  due 
to  the  magnitude  of  the  potential  chosen  for  the  step  response 
(—0.05  to  —0. 10  V).  Dynamic  simulations  for  a  step  in  potential 
from  —0.25  to  —0.30  V  were  also  performed  using  the  model. 
The  qualitative  behavior  of  the  transients  exhibited  by  the  poten¬ 
tial  and  concentration  inside  agglomerate  and  ionomer  film  was 
found  to  be  same.  But  the  variation  of  normalised  concentration 
of  dissolved  oxygen  was  found  to  be  significant. 

The  time  response  of  agglomerate  current  (/agg)  for  step  in 
surface  potential  from  —0.05  to  —0.10  V  is  shown  in  Fig.  14(a). 
It  can  be  observed  that  agglomerate  current  shows  an  instan¬ 
taneous  jump  and  then  decreases  to  a  steady  value.  That  is,  it 
shows  an  inverse  response.  We  abuse  standard  terminology  a 
little  bit  in  terming  this  phenomenon  as  inverse  response.  In 
standard  control  literature,  inverse  response  is  a  phenomenon 
where  the  initial  derivative  sign  is  different  from  the  sign  of  the 
final  steady-state  gain.  However,  the  phenomenon  we  observe  is 
a  sudden  jump  and  a  subsequent  drop  in  the  value,  albeit  with 
the  same  sign  for  the  initial  derivative  and  the  final  steady-state 
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Fig.  12.  Dynamic  simulation  for  normalized  concentration  of  dissolved  O2  inside  agglomerate,  :  (a)  simulation  from  time  t  =  0  to  1.4  x  10  6  s  and  (b) 
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Fig.  13.  Dynamic  simulation  for  normalized  concentration  of  dissolved  O2  inside  ionomer  film,  film:  (a)  simulation  from  time  t  =  0  to  1.4  x  10  6  s  and  (b) 
simulation  from  t  =  0  to  0.035  s. 


gain.  This  can  be  explained  by  observing  the  plots  in  Figs.  10(b) 
and  12(b).  When  a  step  is  introduced  in  the  surface  potential,  the 
agglomerate  potential  (0agg)  shows  an  instantaneous  response, 
see  Fig.  10(b).  However,  the  dissolved  concentration  of  oxygen 
does  not  show  instantaneous  response.  Hence,  a  higher  agglom¬ 
erate  potential  and  high  concentration  of  dissolved  oxygen  are 
the  reason  for  instantaneous  jump  in  the  agglomerate  current. 
Subsequently,  the  dissolved  concentration  of  oxygen  decreases 
to  the  new  steady-state  values  corresponding  to  0r  =  — 0.10  V, 
see  Fig.  12(b).  This  reduced  concentration  and  time  delay  causes 
the  agglomerate  current  to  settle  down  to  its  new  steady- state 
value.  In  Fig.  14(a),  the  magnitude  of  the  inverse  response  is 
small  due  to  the  magnitude  of  the  surface  potential.  A  similar 
inverse  response  in  agglomerate  current  is  observed  for  step  in 
surface  potential  from  —0.25  to  —0.30  V,  Fig.  14(b).  The  effect 
of  potential  magnitude  on  the  magnitude  of  inverse  response 
can  be  clearly  seen  in  this  plot.  In  order  to  check  whether  the 
simulation  results  are  correct,  the  steady- state  values  of  the  ag¬ 
glomerate  current  in  Fig.  14  are  compared  with  those  predicted 
in  the  i-v  curves,  Fig.  5.  The  values  are  found  to  be  same. 

The  dynamics  of  agglomerate  current  in  response  to  a  step  in 
partial  pressure  of  oxygen  at  the  surface  are  illustrated  in  Fig. 


15.  The  dynamics  are  studied  at  two  different  surface  potentials, 
—0.15  and  —0.30  V.  For  both  cases,  simulation  results  show  that 
the  agglomerate  current  increases  to  a  new  steady-state  value 
with  some  dynamics.  In  addition,  the  plots  show  that  at  higher 
potential,  the  dynamics  of  the  agglomerate  current  is  faster,  i.e., 
it  reaches  steady-state  in  lesser  time. 

The  dynamic  results  corresponding  to  step  in  partial  pres¬ 
sure  of  water  vapor  at  the  surface  (/?H?o,r)  are  given  in 
Fig.  16.  When  the  surface  partial  pressure  of  water  vapor 
is  increased,  agglomerate  current  decreases  to  a  new  steady- 
state  value  showing  inverse  response.  This  inverse  response 
can  be  explained  using  the  reason  provided  for  the  steady- 
state  result  in  Fig.  6(b).  When  a  step  is  introduced  in  the 
partial  pressure  of  water  vapor  at  the  surface,  the  water  con¬ 
tent  of  the  ionomer  increases.  This  in  turn,  results  in  a  higher 
ionomer  conductivity  (/cmem)-  Due  to  high  conductivity,  the 
potential  drop  from  the  ionomer  surface  to  the  agglomerate- 
ionomer  interface  decreases.  This  change  in  potential  is  in¬ 
stantaneous,  which  is  exhibited  as  a  sudden  drop  in  the  cur¬ 
rent,  see  Fig.  16(a)  and  (b).  However,  the  dissolved  concen¬ 
tration  of  oxygen  increases  slightly  to  a  new  steady  value, 
but  with  a  time  delay.  This  causes  the  agglomerate  current 
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Fig.  14.  Transients  in  agglomerate  current,  /agg,  for  step  in  </>r:  (a)  step  in  0r:  —0.05  to  — 0.10  V  and  (b)  step  in  0r:  —0.25  to  —0.30  V. 


122 


R.  Madhusudana  Rao,  R.  Ren g asw amy  /  Journal  of  Power  Sources  158  (2006)  110-123 


< 

2.88  r 

2.86 

O) 

2.84 

CD 
_ TO 

r* 

■*-> 

2.82 

C 

<D 

i_ 

2.8 

k. 

3 

O 

2.78 

<U 

•*-> 

CO 

2.76 

u_ 

<D 

E 

2.74  - 

0 

O) 

2.72 

U) 

< 

2.7 

(a) 

2.68  E 
0 

10 


-10 


<)>  =- 0.15  V 

Step  in  pQ  :  0.40  to  0.50  atm 

2 

0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04 

Time  (s) 


(b)  Time  (s) 


Fig.  15.  Transients  in  agglomerate  current,  /agg,  for  step  in  po^y  (a)  step  in  po2y  0.30-0.40  atm  and  (b)  step  in  po^y  0.30-0.40  atm 


to  increase  to  its  new  steady-state  value  with  the  same  time 
delay. 

From  the  y-axis  in  Fig.  16(a)  and  (b),  it  can  be  observed  that 
the  variation  in  current  is  small.  This  corroborates  the  plot  shown 
in  Fig.  6(b),  which  shows  that  the  partial  pressure  of  water  vapor 
has  negligible  effect  on  agglomerate  current. 

The  transients  shown  for  agglomerate  current  (/agg)  for  step 
changes  in  surface  potential  are  typical  and  are  also  observed 
in  dynamic  response  of  fuel  cells.  This  kind  of  behavior  is  also 
shown  by  Wang  and  Wang  [25].  They  have  shown  the  dynamic 
response  of  average  cell  current  density  to  the  step  changes  of 
cell  voltages.  The  plots  show  surge  in  the  current  densities  that 
occur  in  time  scales  that  are  of  the  same  order  of  magnitude  as 
shown  in  Fig.  14.  Our  ultimate  goal  is  to  use  these  dynamic  mod¬ 
els  to  study  operational  problems  in  PEM  fuel  cells  and  stacks. 
Significant  work  has  been  done  earlier  in  our  research  group  for 
phosphoric  acid  fuel  cells  (PAFC)  using  the  dynamic  models  of 
spherical  agglomerates.  It  was  shown  that  these  models  can  be 
used  for  identification  of  problems  during  the  operation  of  large 
sized  PAFC  cells  [26,27]. 


All  the  simulation  results  presented  above  using  steady- state 
and  dynamic  model  of  spherical  agglomerate  can  be  extrapolated 
to  understand  the  behavior  of  PEM  fuel  cell  cathode.  Since  it 
was  convenient  to  carry  out  the  numerical  simulations,  oxygen 
electrode  was  assumed  as  the  reference  electrode.  The  same 
studies  can  also  be  done  considering  that  the  potential  of  the 
solid  phase  inside  the  agglomerate  is  measured  against  a  standard 
hydrogen  electrode  (SHE).  In  such  a  case,  the  local  overpotential 
inside  the  agglomerate  (77 r)  is  defined  as 

T]v  —  Eagg  —  !4gg,e  —  0agg  (24) 

where  Vagg  is  the  potential  of  the  solid  phase  inside  the  agglomer¬ 
ate  measured  against  SHE.  The  subscript  e  denotes  equilibrium 
conditions.  The  i-v  characteristic  curve  and  other  trends  shown 
in  the  plots  for  parametric  studies  will  not  change  if  simulations 
were  carried  out  using  Eq.  (24)  instead  of  Eq.  (13b).  Moreover, 
all  the  trends  that  are  shown  in  the  plots  for  dynamic  studies 
will  also  remain  the  same.  Hence,  all  the  above  results  can  be 
extrapolated  and  used  as  a  basis  to  understand  the  steady- state 
and  dynamic  behavior  of  fuel  cell  cathode. 


Fig.  16.  Transients  in  agglomerate  current,  7agg,  for  a  step  in  pn2o,r:  (a)  step  in  PibOj-  0.10-0.20  atm  and  (b)  step  in  pn2oy  0.10-0.20  atm. 
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4.  Conclusions 

A  dynamic  model  for  spherical  flooded-agglomerate  has  been 
described  in  this  work.  The  flooded-agglomerate  consists  of  a 
thin  film  of  ionomer  covering  a  uniform  mixture  of  carbon  sup¬ 
ported  platinum  and  ionomer.  The  model  takes  into  considera¬ 
tion  the  conservation  of  dissolved  oxygen  and  hydrogen  ions. 
The  complete  characteristic  curve  of  the  agglomerate  is  pre¬ 
dicted  by  the  model.  The  influence  of  surface  boundary  con¬ 
ditions  and  model  parameters  on  the  i-v  curve  were  presented. 
In  addition,  dynamics  of  agglomerate  potential,  dissolved  con¬ 
centration  of  oxygen  and  agglomerate  current  for  step  input  in 
surface  boundary  conditions  were  also  presented. 

In  summary,  the  significant  findings  of  this  work  are: 

•  potential  variation  inside  the  flooded  agglomerate  and 
ionomer  film  is  important  as  significant  gradients  exist  in¬ 
side  the  agglomerate  and  ionomer  film, 

•  with  the  help  of  dynamic  model  studies,  it  was  shown  that  the 
time  required  to  reach  steady  state  for  dynamics  of  dissolved 
oxygen  concentration  and  potential  differ  by  several  orders 
of  magnitude  (10-2  and  10-6  s,  respectively).  The  difference 
in  time  scale  in  the  transients  of  potential  and  concentration 
seems  to  suggest  that  the  dynamics  of  agglomerate  potential 
can  be  neglected  and  the  overall  dynamics  of  the  agglomerate 
model  is  governed  by  the  dynamics  of  dissolved  concentra¬ 
tion  of  oxygen,  co2,agg  and  co2,fiim 

•  the  dynamic  response  of  agglomerate  current  (/agg)  to  step 
changes  in  surface  boundary  conditions  have  been  high¬ 
lighted.  The  time  required  to  reach  steady  state  for  the  ag¬ 
glomerate  current  is  strongly  influenced  by  the  potential  at 
the  ionomer  film  surface. 

All  the  above  results  and  effects  have  been  studied  and  de¬ 
scribed  at  a  single  agglomerate  level.  And  moreover,  the  influ¬ 
ence  of  boundary  conditions  and  parametric  studies  have  been 
done  one  at  a  time,  i.e.,  without  considering  the  effect  of  one 
boundary  condition  on  another  and  one  model  parameter  on  an¬ 
other.  In  the  case  of  complete  PEM  fuel  cell,  this  would  not  be 
the  scenario.  There  would  be  averaging  effects  of  agglomerates 
within  the  reaction  layer,  influence  of  gas  phase  transport  and 
liquid  water,  and  boundary  conditions  that  cannot  be  treated  in¬ 
dependently.  These  effects  are  currently  being  investigated  by 
our  group  as  a  part  of  the  work  on  dynamic  model  for  PEM  fuel 
cell  cathode. 
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